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A model system inspired by recent experiments on the dynamics of a folded protein under the 
influence of a sinusoidal force EM! is investigated and found to replicate many of the response 
characteristics of such a system. The essence of the model is a strongly over-damped oscillator 
described by a harmonic restoring force for small displacements that reversibly yields to stress un¬ 
der sufficiently large displacement. This simple dynamical system also reveals unexpectedly rich 
behavior—exhibiting a series of dynamical transitions and analogies with equilibrium thermody¬ 
namic phase transitions. The effects of noise and of inertia are briefly considered and described. 
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I. INTRODUCTION 

Investigations by Zocchi and collaborators hhh of the 
dynamics of folded proteins under the influence of a si¬ 
nusoidally modulated force provide insight into the me¬ 
chanical response of those molecules to external forces. 
These investigations probe both the interactions and the 
dynamics associated with conformational changes in pro¬ 
teins. The latter aspect of the work promises to enhance 
our understanding of the action of proteins, since con¬ 
formational adjustments, particularly substantial alter¬ 
ations in the tertiary structure of these molecules, are 
central to their action in key biological settings 0111 - 

In the experimental system studied by Zocchi et. al. 
a collection of guanylate kinase enzymes tether 20 nm 
gold nano-particles to a planar gold substrate through 
Cysteins that are mutagenically introduced into those en¬ 
zymes. An oscillating electrophoretic force, generated by 
applying an AC voltage between the gold substrate and a 
parallel electrode, drives the charged gold nano-particles; 
the motion of those nano-particles is then measured with 
the use of evanescent wave scattering. The large number 
of gold nano particles in the sample that participate in 
the motion allow for the detection of displacements that 
are considerably smaller than thermal motion would pre¬ 
dict. Those displacements map directly onto the defor¬ 
mation of the enzymes. 

A key finding that emerges from these investigations is 
the existence of a relatively abrupt crossover as the driv¬ 
ing force increases, from elastic response to a response 
that is dominantly viscous. This crossover is described as 
a “viscoelastic transition.” Reference 0 discusses the na¬ 
ture of the transition and describes two response regimes. 
First, when the displacement is small, the most accu¬ 
rate underlying model of the enzyme appears to be a 


Hookean spring, for which the equation of motion is sim¬ 
ply x(t) = f{t)/k , where f(t) is the applied sinusoidal 
force. This representation of the protein characterizes 
its deformation in terms of a single collective coordinate, 
x(t), and neglects both inertial and dissipative effects. 
On the other hand, for sufficiently large driving force, 
the motion is described in terms of a Maxwellian model 
of a dissipative system [7j, schematically displayed in Fig. 
|T| The equations of motion that govern this system are 

Maxwell 





FIG. 1. Two simple models for a viscoelastic system, along 
with a model for a system governed entirely by viscosity. The 
springs have spring constants equal to k, and the viscous co¬ 
efficient for the dash pots is 7 . 
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where x(t) once again tracks the deformation of the en¬ 
zyme, which is directly related to the displacement of the 
gold nano particle, while y(t) is an “internal” degree of 
freedom of the system. (A further time differentiation 
of Eq. (1.2) renders it in the familiar form relating strain 
rate to the stress and its rate of change.) According to 
[3], the mechanism that drives the transition is an un¬ 
derlying energy versus displacement curve that changes 
from quadratic—i.e., harmonic—at low displacement, to 
linear at higher displacements. The precise nature of this 
form of potential energy curve is not expanded on in that 
reference. 


Inspired by the above results, we have explored the dy¬ 
namical response of a somewhat different— but unified— 
depiction of the driven protein system. Our approach 
is based on a Kelvin-Voigt model for a viscoelastic sub¬ 
stance, as shown in Fig. |T] The spring in this system 
is not, however, strictly Hookean. Rather the energy of 
the spring as a function of displacement, V(x), which 
gives rise to a restoring force F(x), is schematically dis¬ 
played in Fig. [5] In our initial analysis, we ignore inertial 
effects, under the assumption that the system is heav¬ 
ily over-damped. We find that the system displays an 
unexpectedly rich range of behavior, including symme¬ 
try breaking—and restoration—dynamical phase transi¬ 
tions, as well as noise driven rounding and “switching” 
in bi-stability and many features reminiscent of equilib¬ 
rium thermodynamic phase transitions, such as spinodals 
and multicriticality. Preliminary investigations of the ef¬ 
fects of inertia on the dynamical equations reveal an even 
richer range. 

There is at least one precedent for a dynamical tran¬ 
sition in a driven system of the type that is explored 
here. The Suzuki-Kubo equation describes the behavior 
of a mean field version of the Ising model with dissipa¬ 
tive dynamics in which the spin variables are driven by 
a magnetic field with sinusoidal time dependence [8j. At 
low temperatures, the response to the driving field under¬ 
goes a dynamical transition, from an oscillation about an 
equilibrium ferromagnetic state at small amplitude of the 
drive to an oscillation at larger drive amplitudes centered 
about spin magnitude equal to zero [9). This dynamical 
transition can be either continuous or first order, depend¬ 
ing on the temperature. However, the physics underlying 
that model differs fundamentally from the phenomena 
explored here. 

The remainder of this paper is organized as follows. In 
Sec. HI! the basic model is introduced and the determin¬ 
istic, noise free, dynamical phase diagram is displayed 
summarizing the basic behavior. In Sec. m to provide 
additional insight, we provide some analytical results for 
the model with a piecewise continuous restoring force. In 
the following section, Sec. m we make contact with typ¬ 
ical dynamical responses, as measured in experiments on 
viscoelastic materials in general and in the experiments 
of Zocchi et al. in particular. For the basic over-damped 
deterministic case we discuss the nature of the transitions 
and compare the analysis to that of standard mean field 


theory for thermodynamic phase transitions in Sec. [V] 
Section |VI| contains a preliminary investigation into the 
effects of inertia, while the effect of noise on the transi¬ 
tions and response functions is discussed in Sec. m via 
Langevin over-damped dynamics and a master equation. 
Concluding remarks follow, including a short discussion 
of the consequences of a restoring force without the sym¬ 
metry shown in Fig. [2] Appendices contain some details 
on our “standard model,” a comparison to mean field 
thermodynamics and a brief commentary on an alternate 
version of the restoring force. 


II. MODEL, DYNAMICAL PHASE DIAGRAM, 
AND CHARACTERIZATION OF RESPONSE 


Our model is a highly simplified depiction of a folded 
protein, in which we single out a collective degree of 
freedom that we assume dominates the response of the 
molecule to an external force. As a caricature of its much 
more complex structure, we assume a system described 
by the Kelvin-Voigt model of viscoelasticty ej ns , as 
shown in Fig. |TJ which, as we will see, appears natu¬ 
ral for a driven, over-damped, nonlinear oscillator. The 
spring in this model can be thought of as a gross simpli¬ 
fication of the Tirion model for protein interactions El, 
with the caveat that our coordinate, x, is a collective one, 
while the Tirion model replaces the various interactions 
between the actual constituents of a protein by harmonic 
bonds between idealized point-like entities. In addition, 
we assume the possibility of the “cracking” of this pro¬ 
tein under sufficient external stress [121 - 1741 through the 
reversible detachment of the spring. This means that 
the harmonic potential energy expression holds only in 
a limited range of values of x. According to the model 
we adopt, at sufficiently large values of |x| the poten¬ 
tial energy and associated restoring force depart from 
strict Hookean form to approach either a constant en¬ 
ergy (and, correspondingly, no restoring force) or a fixed, 
constant restoring force corresponding to a linear, rather 
than quadratic, energy. Figure [2] is a schematic depic¬ 
tion of a particular version of such a force and associated 
energy function. Inside a region centered at the origin, 
x = 0, the potential is nearly harmonic and the force 
is approximately Hookean. Outside this central region 
the force approaches constant values ±F 0 , and the asso¬ 
ciated confining potential is linear. The precise form of 
the restoring force utilized in our calculations is described 
in Appendix [A] 

We assume initially that the motion of our system is 
highly over-damped, so that the time rate of change of 
the displacement coordinate x is directly proportional to 
the force generating change in the system. Given the 
restoring force shown in Fig. [2] and a periodic external 
force, the equation of motion takes the form 

^ = -[F(x(t)) + AsinM)] (2.1) 

at 7 
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FIG. 2. The potential energy, V(x), (red, dashed) and associ¬ 
ated restoring force, F(x), (blue, solid) that we assume in our 
model. The force approaches ±Fo for large x, as indicated in 
the figure. 


where the parameter 7 encodes viscous effects. 

Figure [3] indicates the resulting behavior when the 
drive amplitude, A , and Fo, the asymptotic absolute 
strength of the confining force, are scanned. In Regions 


1/A 



FIG. 3. Dynamical phase diagram of the types of response to 
a sinusoidal drive exhibited by the equation of motion ( 2.11 
with restoring force F{x) as shown in Fig. [ 2 ] The horizontal 
axis is the absolute value of the asymptotes of the restoring 
force at large \x\. The vertical axis is the inverse of the drive 
amplitude, A. The parameter a controlling the transition 
from Hooke’s law to constant restoring force (see Eqs. (A1 1 
(A 6 l and Fig. 291 has been set equal to 0.01. The frequency 
of the drive, uj, is equal to 1 / 2 , and the viscosity parameter, 
7, is set equal to one. 


1 and la, color-coded blue, the steady state periodic so¬ 
lution, is symmetric about the origin. In Region 2, color- 
coded red, the dynamically stable solution is skewed, ei¬ 
ther to the right or the left of the origin. As it turns 
out, skewed solutions appear in pairs; see Appendix [B] 
There is also a symmetric steady state solution that is, 
however, dynamically unstable. In Region 3, color-coded 
green, three dynamically stable, steady state solutions 
exist, one symmetric and the other two skewed. In addi¬ 
tion, there are two skewed, dynamically unstable, steady 
state solutions. Finally, a region complementing Region 
3, separating Region 1 from Region 2, is not shown as it 
is exceedingly narrow and beyond the resolution in Fig. 

El 


Figures |4]-[6] illustrate the (steady state) solutions char¬ 
acteristic of the four regions in Fig. [3j The time interval 
shown corresponds to two periods of the driving force. 
The parameters A and Fo corresponding to the region 
of the phase diagram are provided in the captions, while 
k = 2 , corresponding to the spring constant in the har¬ 
monic regime, and to = 0.5 for the driving force are kept 
constant in these figures. The friction constant 7=1 
sets the time scale. Given the very different properties 


x(t) 



FIG. 4. The steady state response in Region 1 (solid curve) 
and in Region la (dashed curve). Both are dynamically sta¬ 
ble. The inverse of the drive amplitude, 1/A, is 0.53 for the 
Region 1 response and 0.51 for the Region la response. The 
asymptotes of the restoring force outside the harmonic regime 
are ±Fo with Fo = 0.6, and the value of the spring constant 
in the Hooke’s law region is k = 2. The frequency of the 
driving force is ui = 0.5. The transition between a harmonic 
restoring force and a force equal to ±Fo occurs at \x\ ~ 1. 
The restoring force is given in Appendix [A] with parameter 
a = 0.01. 


x{t) 



FIG. 5. The three steady state responses in Region 2. The 
solid, skewed, curves are dynamically stable and the dashed, 
symmetric, curve is dynamically unstable. The quantity Fo 
is 0.1 and 1/A = 0.45. All other parameters are the same as 
in Fig. 0 


of the steady states in the four regions, the boundaries 
between the regions are necessarily sharp. We will return 















4 


x{t) 



FIG. 6. The five steady state responses in Region 3. The solid 
symmetric curve and the two solid skewed curves are dynam¬ 
ically stable, while the two skewed curves shown dashed are 
dynamically unstable. The quantity Fa is 0.1, and the inverse 
of the drive amplitude, 1/A, is 0.35. All other parameters are 
the same as in Figs. [4] and [5] 


to the nature of the transitions that take place as those 
boundaries are traversed. 

For the time being, it is worthwhile to consider the two 
solutions for x(t ) shown in Fig. [4j They differ in a few re¬ 
spects. First, although the driving force strengths differ 
by less than 5%, the displacement amplitudes are quite 
different. Second, the curve describing the solution in Re¬ 
gion 1 is nearly sinusoidal (in fact, nearly in phase with 
the drive, which goes as sin(wf), while the steady state 
solution for Region la is somewhat distorted. Further¬ 
more, the latter solution is well out of phase with respect 
to the drive. As we will see, the steady state solution in 
Region 1 is close to elastic—i.e., nondissipative—while 
the the behavior in Region la is strongly dissipative. 


III. ANALYTICAL SOLUTION IN A LIMITING 
CASE 


To gai n a perspective on the nature of solutions of 
Eq. (2.1), we consider the limiting assumption of a sharp 


break between the regime in which the restoring force is 
strictly linear and the range of the displacement variable 
x in which it is considerably gentler. If the restoring force 
outside the harmonic regime vanishes , the function F{x) 
has the form 


F{x) = 


—kx |ir| < Xo 
0 lad > Xo 


(3.1) 


The solutions to Eq. (2.1) in the two regimes are 

Xl {t) = Be~ kth 

A ( k \ 

- . — cos (tut + arctan — ) (3.2) 

^J 2 UJ 2 + k 2 V 7W/ 


x e (t) = B' -cos(wf). 

7 ui 


(3.3) 


where the subscripts i and e refer to behavior in the “in¬ 
terior” region \x\ < Xq and the “exterior” region |x| > Xq- 
A complete solution to the equation of motion requires 
adjusting the coefficients B and B' so as to match 27 (f) 
and x e (t) at the boundary between the two regimes. 

When the frequency of the driving force is small, the 
characters of the “exterior” and “interior” responses dif¬ 
fer fundamentally. In the exterior region, the steady state 
velocity, v e (t) = dx e (t)/dt, is in phase with the driving 
force, which means that the drive’s energy feeds opti¬ 
mally into viscous damping. On the other hand, when 
70 j -C k, dxi(t)/dt is almost ninety degrees out of phase 
with the driving force, so the response is nearly the same 
as at static equilibrium, in that the two terms in square 
brackets in ( 2 . 1 ) nearly cancel, and the dissipation is 
quite small. In the regime of low frequency response, 
we can roughly characterize the dynamics in the exterior 
region as viscous and the behavior in the interior region 
as elastic. It is therefore reasonable to expect that as the 
driving amplitude, A, increases, the response will evolve 
from elastic to viscous. 

However, as strongly indicated by the phase diagram 
shown above, the change from elastic to viscous response 
entails an abrupt transition in the dynamical behavior. 
This transition occurs near the drive amplitude for which 
the coordinate x visits the region outside the harmonic 
regime (in fact, it occurs at a slightly lower drive am¬ 
plitude). At this point skewed solutions to the equation 
of motion appear and, in fact, coexist with solutions re¬ 
maining entirely within the harmonic regime. 

Figure [7] shows two full p erio ds of two skewed solution 
to the equation of motion ( 2 . 1 ) with confining potential 
& The displacement limits outside of which which 


x(t) 



FIG. 7. Two skewed solutions to Eqs. ( |2.1| ) and ( |3.1| ). The 
value of xo has been set equal to one. 


the restoring force vanishes are indicated by dashed lines. 
Note the slope discontinuities when those limits are tra¬ 
versed. This is permissible in a system in which inertia 
is ignored. The trajectories consist of the two solutions 
(3.2) and (3.3), grafted together at the boundaries be¬ 
tween their regimes of applicability. From the nature of 
the solution for |a;(f)| > xq, we see that, in the case of the 
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solution skewed above the x axis, the intervals in which 
the solution for \x(t)\ < Xq applies, as depicted in Fig. [7j 
will satisfy 

(2 n + l) 7 r/w + to<t< (2 n + 3 ) 7 r/w — to , (3.4) 

corresponding to time windows lying symmetrically 
within the interval between successive odd multiples of 
the period of the forcing term, 2'k/lo. We can write for 
the form of a solution in such a time window 


v{t) 


A 


WsJ 1 + (fc/w) 2 

_|_£ig— k(t— (2n+l)7r/w— 1 0 ) 


cos (cot + arctan(fc/w)) 


(3.5) 


Here and henceforth, we have set the parameters 7 = 
l,:ro = 1- hr order for the solution in (3.5) to properly 


match the solution for |:r(f)| > xq at the limits of the 
window, we require 


A 


x 0 = - 


x 0 = - 


cos(( 2 n + l) 7 r + wfo 


Wi/l + (fc/w ) 2 

arctan(fc/w)) + C 

- cos((2 n + 3)tt - ujt 0 

wa /1 + {k/uiy 

arctan(fc/w)) + Ce~ 27lk/uj+2kto 


(3.6) 


(3.7) 


We can eliminate the constant C between these two equa¬ 
tions, and we are left with the equation for to 


A 


cos (—Luto + arctan(fc/w)) 


UJy/l + ( k/0J ) 2 
_ e ~ 2 Trk/ui+ 2 kt 0 cos ( u t 0 arctan(fc/o;)) 

+x 0 (e- 2 * k '“ +2kta - l) = 0 


(3.8) 


The left hand side of (3.8) as function of t 0 is displayed 
in Fig. |8j with the variables xq, k and ui set so as to 
correspond to their values in the phase diagram in Fig. [3] 
As shown in the figure, smaller amplitudes correspond to 
lower curves. Note that there is always a solution uto = 
7 r, which says that the window in (3.4) has shrunk to zero. 


This is reasonable since it is always possible, by setting 
initial conditions appropriately, to construct a solution 
in which x(t) lies completely in the regime of vanishing 
restoring force. When to = tv/oj is the only solution to 
(3.8), then solutions to the equations of motion (2.1) with 


F(x) given by (3.1) either lie completely in the regime 
| a: | < xq (where, in this example, Xq = 1 ) or entirely 
outside of it. At a threshold value of the amplitude, two 
more apparent solutions arise, as for the solid curve in 
Fig. m 

The next step is to determine whether the new inter¬ 
sections correspond to solutions of the equation of mo¬ 
tion, and, if so, whether those solutions are dynamically 
stable. As it turns out, one of the two new solutions 



to 


FIG. 8. The left hand side of (3.81 as a function of to 
xq = 1, k = 


with 


2 and uj = 0.5. The three curves correspond to 
the following three values of the drive amplitude: A = 1 (blue 
long-dashed curve), A = 2 (red dashed curve), A = 3 (black 
solid curve). 


of (3.8) -the solution corresponding to the larger value 


of to —does indeed correspond to a legitimate response. 
Furthermore, the response is dynamically stable. Ap¬ 
pendix [C] describes the analysis. 

Another interesting aspect of the dynamical behavior 
is the existence of an exceedingly small range of ampli¬ 
tudes, A , in which a dynamically stable, skewed solution 
coexists with a symmetric solution lying entirely inside 
the Hooke’s law region. As an example of coexistence 
Figure |9] shows the left hand side of Eq. (3.8) when 



to 



to 


FIG. 9. The left hand side of Eq. ( |3.8| ) when k = 2, uj = 0.5, 
7 = 1, xo = 1 and the drive amplitude A is equal to y/k 2 + u 2 . 
The lower portion of the figure shows the expression for the 
entire range of to, from 0 to n/oj. The upper portion magnifies 
the portion of the graph in which the left hand side of (3.81 
passes through zero. 


the parameters are adjusted so that there is a solution 
to the equation of motion that just fits in the region 
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—Xq < x < xq. Even though there is such a trajectory, 
an additional skewed and dynamically stable solution to 
the equation of motion with F(x) given by © 

also exists (actu ally a pair, see below), corresponding to 
the solution to ( |3.8| ) in which to ~ 2.654. The region of 
coexistence is spanned as the drive amplitude A changes 
by about one part in 10 6 . Figure 10 shows the three 
stable solutions. Not shown in the figure are the two 


x(t) 



FIG. 10. The three stable solutions to the equation of motion 
(2.11 with restoring force given by (3.11. The parameters are 
the same as in Fig. [9] The two skewed solutions are shown 
as solid curves and the dots indicate the symmetric solution, 
which just fits in the Hooke’s Law region for the restoring 
force. The boundaries of that region are shown as dashed 
lines. 


unstable skewed solutions to the equation of motion that 
lie between the symmetric solution and the stable skewed 
solutions. 


Given a solution, x(t ) of the equation of motion (2.1) 
with the kind of restoring force considered here, i.e. a 
restoring force with the property F(—x) = —F(x), it 
is straightforward to show that the equation of motion is 
also satisfied by —x(t+n/tjo). From this we can infer from 
one skewed solution another one, skewed in the opposite 
direction. This feature was used to construct the pair 
of skewed solutions in Fig |10| Thus, we have confirmed 
the existence of two distinct stable, skewed solutions to 
the equation of motion for sufficiently high values of the 
driving amplitude. See Appendix |B| 


IV. LOSS AND STORAGE COMPLIANCES 


To further quantify the transition, we introduce a gen¬ 
eralization of the way in which the response of a linear 
viscoelastic system is characterized in terms of dissipa¬ 
tive, or loss, compliances and storage, or elastic, com¬ 
pliances. Given the notation we utilize here, the stress 
exerted on the system is Asin(wt). The response of the 


system, if linear, would be of the form 

xit'j ^diss(f') + ^el (f) 

= —X'd cos (cot) + X e sin(wf) 

= — j^Acos^t) + jiAsm(ut ), (4.1) 


where j\ is the storage compliance and j 2 is the loss 
compliance [TO], and the subscripts on the first line of 
(4.1) stand for “dissipative” and “elastic.” 

In the case of nonlinear response, it is possible to gen¬ 
eralize the above decomposition by extracting the contri¬ 
bution to the response that gives rise to dissipation. We 
do this by taking the integral of the steady state solution 
over a period, 


w r 27l/ui 

Yd= — / x(t) cos (ojt)dt , (4.2) 

7T Jo 

in which case we can write 


x(t) = Y d cos(ujt) + (x(t) — Y d cos(u>t)) 

= x diss (t) + x e i(t) (4.3) 

The entirety of the dissipative response is contained in 
Xdi S s(t) as defined above, in that 

f 2n/ul dx At) 

/ sin(wf )—-—dt = 0 (4.4) 

Jo dt 


This means that the quantity —Y d /V2A (see below) plays 
the role of the dissipative compliance, while there is not 
necessarily any quantity that can be unambiguously as¬ 
sociated with the storage compliance. 

We now return to our ’’standard model” shown 
schematically in Fig. [2] Figure |TT] displays the two so¬ 
lutions for x(t) shown in Fig. |4j divided into dissipative 
and elastic components. As is clear from the lower plot in 
the figure, which refers to the response in Region la, the 
elastic response is not simply proportional to the driving 
force, Asin(wt). 

However, based on the developments above, we can 
reduce the elastic and dissipative components of the re¬ 
sponse to two numbers. If 


then 


/>27T / UJ 




2tt 


•:(t) 2 dt , 


(4.5) 


/>27T / cu 

~ / *^diss(^) dt 


UJ 
27 T 


UJ 
27T 


r 2ir / uj 


x e \(t) 2 dt 



(4.6) 


(4.7) 


which means we can define an overall compliance as It /A, 
with a storage compliance, j 1 = I s /A, and a loss compli¬ 
ance, j 2 = Id,/A. 
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•£diss,el (t') 



FIG. 11. The two solutions for x(t ) shown in Fig. [4] broken 
into dissipative contributions (solid curves) and storage con¬ 
tributions (dashed curves); see Eq. ( |4.3| ). The top plot shows 
the solution in Region 1, in which the response is dominantly 
elastic and the bottom plot shows the solution in Region la 
in which the response is more dissipative. Recall that the in¬ 
verse of the drive amplitude is 0.53 in the case of the top plot 
and 0.51 for the bottom plot, and Fq = 0.6. 


As a final technicality, we need to address the case in 
which there are dynamically stable skewed solutions to 
the equation of motion. The measurements that moti¬ 
vated this study are made on a collection of driven oscil¬ 
lators; furthermore, thermal noise is substantial. As we 
will see such a system will equilibriate into an ensemble 
in which the two skewed responses, when they exist, are 
equally represented. This means that under such con¬ 
ditions, we should replace x(t) in the equations defining 
the compliance by 


(x(t)) = (x(t) — x(t + 7r/w))/2 


(4.8) 


Figure [redisplays an example of the case in which the 
response curve is skewed and the average in (4.8) has 
been performed The full response is shown as well as the 
two components that comprise it. 


V. ON THE NATURE OF THE DYNAMICAL 
PHASE TRANSITION 


Our focus has been the steady state behavior of the 


dynamical system described by (2.1) When the crossover 


x{t) 



FIG. 12. The average of the two stable skewed response 
curves. The solid black curve is the full response. The long 
dashed blue curve is the dissipative contribution, going as 
cos(cjt), which we term Xdiss(t). The dashed red curve is what 
remains when the dissipative component is subtracted from 
the full response. We term that response x e \(t) (See (4.31. 
Here, k = 2, lu = 0.5, Fq 0, a = 0.01 and A, the drive 
amplitude, is equal to 1.9. 


from a restoring force described by Hooke’s Law to a con¬ 
stant restoring force F 0 is not abrupt, the methods of Sec. 
|III| cannot be applied. An alternate approach that proves 
powerful, useful and universal, utilizes a standard map 
or recursion relation. This relation follows from the fact 
that the solution of a first order differential equation like 
© is determined entirely by a single initial condition. 
That is, if we know x(0), then we can use the equation 
to forward integrate and determine x(t) at all subsequent 
times. In particular, we can create a map that takes us 
from x(0) to x(2n/ix), the value of the dynamical variable 
one period later. A steady state solution to the equation 
of motion will have the property x(2ir/cj) = a;(0)— if we 
can discount the possibility of solutions with a period at 
a subharmonic of the driving force. Figure |13| d isplays 
such a recursion graph, calculated from Eq. ( |2.1| ) with a 
particular choice of parameter values. Also shown in that 
figure is the 45° line corresponding to x(2n/u>) = a;(0). 
The small open circles indicate intersection of that line 
with the recursion curve, corresponding to steady state 
solutions of the equation. 

We can assess the dynamical stability of the steady 
state solutions in the standard fashion by iterating the 
recursion relation Haul]. The graphical version of this 
process is displayed in Fig. [Mj in which we focus our 
attention on the section of the curve containing the in¬ 
tersections. The stepwise curves between the map and 
the 45° line are graphical renditions of the result of it¬ 
eration of the recursion relation. As indicated by the 
arrowheads on those curves, repeated calculations of the 
quantity x(t) at succeeding intervals of one period tend 
away from the central intersection and towards one of 
the flanking ones. We are led to the conclusion that the 
central solution for a steady state is dynamically unsta¬ 
ble, while the two solutions that flank it are, by contrast, 














FIG. 13. Example of the map co nnec ting a;(0) with x{2ty/lj), 
the quantity x(t) satisfying Eq. (2.11. Also shown is the line 
x(2ty/uj) = *(0). In this case, the drive amplitude, A is 2.75, 
the frequency, u>, of the drive is 2, k is 2, the parameter a 
is 0.15 and the asymptotic force, Fo, is equal to zero. As 
always, xq = 1 and 7 = 1. The small open circles denote the 
intersections between the map and the 45° line corresponding 
to steady state solutions of the equation of motion. 


FIG. 15. The total compliance, It/A, plotted agains the drive 
amplitude, A. Here, Fo=0, ui = 0.5, k — 2 and the parameter 
a = 0.01. Here and in Fig. [l6]the compliances are as defined 
in Sec. |!V] 


HA 



FIG. 14. The result of graphical iterations of the recursion 
relation embodied in the curve displayed in Fig. [13] The 
broken curves impinge on the recursion curve vertically and 
on the 45° line horizontally. 


FIG. 16. The two contributions to the total compliance 
graphed in Fig. [l5] the dissipative compliance Id/A (long 
dashed blue curve), the storage compliance I g /A (dashed red 
curve) and the total compliance, (thin black curve). 


dynamically stable. Such a recursion relation yields the 
solutions one finds in Region 2 of the dynamical phase 
diagram as shown in Fig. [3] 

Given a steady state solution for a particular a;(0), we 
can then compute x{t) for the entire interval between 
t = 0 and t = 27r/w, and from this the response properties 
of that solution. As an example, we can explore in more 
detail the transitions between the various regions shown 
in the phase diagram in Fig. [3] We start by looking at 
the portion of the phase diagram that is exactly on the 
vertical axis, i.e. for which the large-cc asymptote of the 
restoring force is Fo =0. Making use of the definitions in 
Sec. [TV] for total, storage and loss compliances, we obtain 
results for the quantities It/A, I s /A and Id/A. Figure [l5| 
shows the total compliance, It/A, plotted in terms of the 
amplitude, A, of the drive. Note the discontinuity in the 
compliance, which can be taken as evidence for a first 
order transition. Figure [16] shows the two contributions 
to the total compliance, obtained from Id/A, and I s /A. 


Close investigation reveals a very narrow regime 
around the point of discontinuity in which there is coex¬ 
istence between the responses corresponding to the com¬ 
pliances to the right and left of the transition point in 
Figs. [T5]and[T6| The regime is sufficiently thin that it is 
undetectable given the resolution of the figure. 

The discontinuities in the dynamic response graphed 
in Figs. [15] and [16] imply a dynamical transition with 
the characteristics of a first order thermodynamic phase 
transition, in that there are discontinuities in key quan¬ 
tities. However, a transition with discontinuities is not 
inevitable. If the change in the restoring force from har¬ 
monic to a constant (in the example being discussed: 
Fo = 0) is sufficiently gentle, then the dynamical trans¬ 
formation becomes continuous in the thermodynamic 
sense. That is, it remains sharp, but physical properties, 
such as dynamical moduli, no longer exhibit discontinu¬ 
ities in their dependence on the drive amplitude. 
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The change in the character of the transition as the 
properties of the restoring force are altered is illustrated 
in Figs [TT] and [l8j which display plots highlighting the 
intersections of recursion curves with the 45° line. In 
particular, the figures show how the points of intersec¬ 
tion change with drive amplitude. The plot in Fig. [17] 
contains four recursion curves, corresponding to four dif¬ 
ferent drive amplitudes. In this case there is a first order 
dynamical phase transition. The parameters are—with 
the exception of uj —the same as in the phase diagram in 
Fig. m with Fq = 0. The drive frequency has been set as 
uj = 2 in order to make certain features of the curve more 
visible and to expand the coexistence region between Re¬ 
gions 1 and 2. The longer dashed blue curve intersects 
the 45° line once, corresponding to a single symmetric 
and stable dynamical steady state. This curve is charac¬ 
teristic of low drive amplitude and (more nearly) elastic 
response. The (short-) dashed red curve corresponds to 
the drive amplitude at which one sees the onset of four 
more steady state solutions to the equation of motion, 
all of them skewed, two stable and two unstable. In the 
regime represented by the solid black curve, those addi¬ 
tional solutions are clearly visible, and three dynamically 
stable steady state solutions to the equation of motion co¬ 
exist - one symmetric and two skewed. The long-dashed 
red curve corresponds to the high amplitude regime, Re¬ 
gion 2 in the phase diagram, in which there are two stable 
skewed steady state solutions and one unstable symmet¬ 
ric solution to the equation of motion. The restoring force 
is as given in Appendix [A] with parameter a = 0.01, the 
same value as was used to generate the phase diagram in 
Fig. § 



FIG. 17. Illustrating a first order transition between Region 
1 and Region 2. Referring to Appendix [A] the parameters in 
the restoring force and the equation of motion are as follows: 
oo = 2, k = 2, xo = 1, a = 0.01 and Fo = 0. The longer dashed 
blue curve and the shorter dashed red curve correspond to 
Region 1 in the phase diagram in Fig. [3] The longest dashed 
red curve corresponds to Region 2, and the solid black curve 
to a coexistence region between Regions 1 and 2. Such a 
region exists but is exceedingly narrow for the set of variables 
used to generate Fig. [3] 


Figure [18] displays the effect of smoothing out the tran¬ 


sition region in F(x). Here, the parameter a has been 
set equal to 0.15. In this case, the transition from a sin¬ 
gle, stable and symmetric, steady state solution to three 
steady state solutions—two skewed and stable and one 
symmetric and unstable—occurs continuously with no 
coexistence region and no discontinuities in dynamical 
response. The drive frequency, uj, has again been set 
equal to 2, while the spring constant k in the Hooke’s 
Law region is maintained at 2. The short-dashed blue 
curve corresponds to Region 1, the long-dashed red curve 
is characteristic of Region 2, and the solid black curve il¬ 
lustrates the onset of the transition between one symmet¬ 
ric, stable solution to the equation of motion and three 
solutions, two skewed and stable, one symmetric and un¬ 
stable. No region of coexistence separates the two. 



FIG. 18. Illustrating a continuous transition between Regions 
1 and 2 in Fig. [3] when a is increased from 0.01 to 0.15 so as 
to stretch out the transition region between a Hooke’s Law 
restoring force and a vanishing restoring force, i.e. Fo = 0. 
Otherwise the parameters utilized are the same as in Fig. 
|18| See Appendix A. The dashed blue curve corresponds to 
Region 1 and the long dashed red curve to Region 2. The 
solid black curve corresponds to the onset of the transition 
between the two regions. 


Figure [T9] illustrates the compliances associated with 
Fig [18| quantified in terms of those defined in Sec. m 
Because of the particular relative values of uj and fc, the 
compliance in Region 1 of the phase diagram divides 
equally into viscous and elastic, while in Region 2 vis¬ 
cous response increasingly dominates. 

The compliance curves are continuous and without ev¬ 
ident features marking the transition from Region 1 to 
Region 2. Very close inspection reveals singularities in 
the form of mild slope discontinuities in It/A, I s /A and 
Id/A , which are unlikely to be detectable in any exper¬ 
imental realization of this system. Figure 20 shows the 
total compliance, It/A, in the immediate vicinity of the 
transition. 

The recursion curves bear a close relationship to the 
free energy extremum curves that arise in a simple order 
parameter-based model of a tricritical system. To be 
specific, the intersection of the recursion curves with the 
45° line correspond with the intersection of the curves 
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FIG. 19. In the case of parameter values leading to the plot in 
Fig. |18| the total compliance, It./A, (solid black curve) and 
the three contributions to that compliance: the dissipative 
compliance Id/A (long dashed blue curve) and the storage 
compliance I a /A (dashed red curve). Here, the transition 
is continuous. Also shown on the graph is the value of the 
amplitude, A at which the transition between Section 1 and 
Section 2 in the phase diagram occurs. For definitions of the 
compliances see Sec. m 


iti a 



FIG. 20. The total compliance, It/A, in the immediate vicin¬ 
ity of the continuous transition from Region 1 to Region 3 in 
the dynamical phase diagram. 


for dp(ip)/dip with the abscissa, as shown in Figs. [32] 
and [34] in Appendix [D] which reviews the classical order 
parameter based model of a tricritical point. 


VI. ON THE EFFECTS OF INERTIA 


So far we have ignored any effects of inertia on the re¬ 
sponse of the system. In light of the unexpected nature 
of the over-damped response and the typically significant 
effects of inertia, including resonant behavior, it is worth¬ 
while to ask whether the dynamical transitions observed 
in the absence of inertia survive its introduction into the 
equations of motion. To this end, we have performed 


preliminary studies of an extension of ( 2 . 1 ) that incor¬ 


porates an inertial term. The new equation of motion 
is 

^ djr(f) +7 ^1 = F(x(t)) + Asin(u;f), (6.1) 

which can be rewritten as 


dx(t) 

dt 

dpjt) 

dt 


P(t) 

m 

—7 x(t) + + Asin(ujt) 


( 6 . 2 ) 

(6.3) 


The search for steady state solutions to the above equa¬ 
tion of motion can be formulated in terms of a pair of 
recursion relations. That is, on the basis of Eqs. (6.2) 
and (6.3) we can construct a map from Xq = x(0) and 
Po = p(0) to x\ = x(2tt/uj) and p\ = p(2ir/u>). The map 
takes the form 


x 1 =X(x 0l p 0 ) (6.4) 

Pi = P{x 0l p 0 ) (6.5) 


Steady state solutions to the equation of motion are Xf 
and pf, where 

x f= x i x f,Pf ) ( 6 - 6 ) 

Pf = p { x f,Pf) ( 6 - 7 ) 


Taken separately each of the relations above generates 
a curve in the ( 2 /, pf ) plane. Figure 21 displays two 


such curves. The point at which the two curves intersect 
corresponds to a simultaneous solution of ( 6 . 6 ) and (6.7) 
and thus a steady state solution of (6.2) and (6.3). The 
solution x[t) arising from that intersection, which corre¬ 
sponds to a relatively weak driving force, is displayed in 
Fig. [22j Also shown in that figure are the nominal limits 
of the harmonic restoring force, at x = ± 1 . 

By contrast, when the amplitude of the drive increases, 
more than one solution to the fixed point equations ( 6 . 6 ) 
and (6.71 appear, as shown in Fig. 23 The curves for x(t) 
corresponding to the three intersections are shown in Fig. 
[24| The relationship between the two skewed solutions 
is the same as the relationship between “mirror image” 
solutions established for solutions to (2.1) in Appendix 
|B| In fact, the argument in Appendix B is easily extended 
to incorporate the inertial term in the equation of motion 

The question of the stability of the solutions of ( 6 . 6 ) 
and (6.7) can be addressed in a straightforward manner. 
Imagine that x f and pf are such simultaneous solutions. 
Then, if xq = Xp + Ax and Po = Pf + A p, 


x\ = xf + Ax' 

= X{x F + A x,p F + A p) 

= xf + X x Ax + X p Ap (6-8) 

Pi = Pf + A p' 

= Pf + P x Ax + PpAp (6-9) 


or 


( Ax' 
\Ap' 


X x X p \ ( Ax \ 
P x P p ) \Ap J 


( 6 . 10 ) 
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Pf 



x f 


FIG. 21. A curve corresponding to ( |6.6| ) (solid curve) and 
(6.71 (dashed curve). In this case, we have taken k = 2, 
7 = 1, a = 0.01 in the expression for the restoring force 
(see Appendix [A| and the force outside the Hookean region, 
Fo = 0. The effective mass is m = 1, and the drive amplitude, 
A = 1.82. 


x(t) 



FIG. 22. The solutions, x(t), corresponding to the point of 
intersection in Fig. [2l] The response is entirely inside the 
regime in which the restoring force is Hookean. 



X f 


FIG. 23. Another pair of curves corresponding to (6.61 (solid 
curve) and (6.7) (dashed c urve) . In this case the parame ters 
in the equation of motion (6.1l are the same as in Fig. |21| 
except for the drive amplitude, A, which is now 1.87. Note 
the presence of three intersections of the two curves. 


x(t) 



FIG. 24. The three solutions corresponding to the intersec¬ 
tions in Fig. |23| The two solid curves are the two stable, 
skewed solutions to the equation of motion (6.11, while the 
dashed curve is a symmetric and unstable solution to that 
equation. The horizontal lines indicate the limit of the region 
in which the restoring force is Hookean. 


from which we infer that the stability of the solutions is 
controlled by the eigenvalues of the matrix on the right 
hand side of (6.10). If both have absolute value less than 
1, the fixed point at Xf, yF is stable; otherwise it is 
unstable. It is possible to determine the elements of the 
matrix numerically by exploring the recursion relations 
in the vicinity of the fixed point. Using this analysis, we 
find that the single intersection shown in Fig. [2l]is stable 
and that the intersections on the far right and far left of 
Fig. [23] also represent dynamically stable solutions to the 
simultaneous equations (6.6) and (6.7) while the central 


intersection, which gives rise to the steady state solution 
shown dashed in Fig. 24, corresponds to a dynamically 
unstable solution to the equation of motion (6.1). 

An alternate approach to the stability analysis of the 
steady state solutions to (6.1) is described in Appendix 


The above results strongly indicate the the types of dy¬ 
namical transitions that we find when inertia is neglected 
are also to be expected in a system that manifests inertia. 
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A more comprehensive exploration of the nature of vis¬ 
coelastic behavior in such systems, including the range of 
dynamical transitions that one may encounter, remains 
to be undertaken. 


VII. NOISE 


A. Langevin equation approach 


As noted in the Introduction our original motivation 
was a series of experiments on the mechanical properties 
of a driven natively folded protein in solution at room 
temperature. Clearly this is a noisy environment, and 
it behooves us to examine which features we have un¬ 
covered survive the introduction of fluctuations, thermal 
and otherwise. In this section we briefly and qualitatively 
consider some of the modifications that arise when noise 
is introduced into the over-damped equation of motion, 
Eq. (2.1l. For a specific example, we take the nonlin¬ 


ear restoring force to be the one defined in Appendix A, 
with Fq = 0, a = 0.01. The restoring force is linear in 
displacement for |cc| less then about unity, decaying to 
zero restoring force outside that interval, as sketched in 
Fig. [2] with F 0 = 0. 

With noise the dynamical equation is 

= ~[~V + Asin(cjt)] + r](t) (7.1) 

where rj(t) is a Gaussian random variable, corresponding 
to white noise, with correlator 


(v(t)ri(t')) = TS(t — t') 


(7.2) 


We do not assume a fluctuation-dissipation relation and 
take 7 and T as independent parameters. 

A simple Euler forward integration leads to 

x(t + St) = x(t) + — [— V'{x{t)) + Asin(wf)] + VTStyi 
7 

( 7 - 3 ) 

where yi is a Gaussian random variable with zero mean 
and unit standard deviation. As an example, we choose 
parameters for which the noise-free system has a first 
order transition from a “phase” at sufficiently low driv¬ 
ing amplitude A in which there is a single, stable steady 
state, which is symmetric, to a phase (probably not phys¬ 
ically reachable, occupying a tiny region of parameter 
space) at higher driving amplitude with three stable so¬ 
lutions, one symmetric and a pair of skewed solutions 
and, finally, at even higher amplitude drive to a phase in 
which only the skewed solutions are stable. For reference, 
for the parameters we choose, the noise-free dynamical 
transition occurs with amplitude about A ~ 1.89. A re¬ 
alization of the three-solution situation with added noise 


is shown for A = 1.89 in Fig. 25 


The noise level T = 1CP 4 has been chosen so that the 
noise-driven switching between the stable solutions is ap¬ 
parent. The nature of the switching will be considered 


x(t) 



FIG. 25. One realization of the solution to the stochastic 
equation (7.1 1 when the noise-free system is in the vicinity of 
a first-order dynamical phase transition from a phase with a 
single symmetric solution to a phase with three stable (noise 
free) steady state solutions, one symmetric about x = 0 and 
a pair of skewed solutions. The restoring force is provided 
in Appendix |A| with a = 0.01,/c = 2,Fq = 0. Other param¬ 
eters: drive amplitude, A = 1.89, frequency, to = 0.5, noise 
correlator, T = 0 . 0001 . 


in future work. With this level of noise, the single sym¬ 
metric solution below the transition is unremarkable and 
not shown. 

For modest levels of noise, such as those considered 
above, qualitatively one expects the sharp noise-free tran¬ 
sitions to be rounded, but to retain hints of the underly¬ 
ing noise-free transitions. This is, in fact realized. How¬ 
ever, the calculation turns out to be more straightfor¬ 
wardly carried out in the context of the Fokker-Planck 
equation, to which we turn now. 


B. Master equation 

The dynamical equation with noise, Eq. |7.1| is Marko¬ 
vian, of the form x = h(x , t) + rjit), where, as above, rj(t) 
is delta-correlated with constant T, and h(x,t) contains 
the “systematic” restoring force and the time-dependent 
driving force. As shown, e.g., in m, one can reformulate 
the analysis of this equation into a probability distribu¬ 
tion, P(x,t), of an ensemble of solutions to the noisy 
equation of motion. This distribution is governed by the 
Fokker-Planck equation, 

P{x,t) = D {1) (x,t)P{x,t ) 

+ (j^) D^(x,t)P(x,t) (7.4) 

D^\x,t) = h(x,t) (7.5) 

D ( - 2 \x,t) = T/2 (7.6) 

In our case, the quantity h(x,t ) is the driving force 

F(x(t)) + Asin(wf). As an example of the utilization 
of this equation, we apply it to the case of a discontin¬ 
uous transition described in Sec. [V] Initially, we take 
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FIG. 26. The probability distribution, P(x, t) as given by 
the solution of the Fokker-Planck equation (7.41 (7.61, with 


A = 2.7, w = 0.5, k = 2, Fo = 0 and T = 0.2. Ten periods 
of the driving force are shown. Lighter colors correspond to 
higher values of P(x,t). The distribution at t = 0 was chosen 
to center on one of the dynamically stable skewed solutions 
to the deterministic equation of motion. 


amplitude there are two stable skewed solutions and one 
unstable symmetric solution. The initial distribution was 
chosen to center on one of the skewed solutions. As il¬ 
lustrated in Fig. [26] the distribution evolves toward one 
that is symmetric. That is, both skewed solutions are 
equally represented in the steady state distribution. In 
fact, one can r eadil y dem onstrate that if P{x , t) is a so¬ 
lution to Eqs. (7.4)—(7.6), then so is P(—x,t + n/uj). In 
Fig. [27] we show the steady state distribution for one 
period, along with the two stable skewed solutions and 
the unstable symmetric solution for zero noise. 


C. Compliances in the presence of noise 


Making use of the steady state probability distribu¬ 
tion, one can calculate the effect of noise on compliances. 
Figure [28] displays the results for the compliances in the 
case in which the parameters in the noise-free equation 
of motion are such that Eq. © predicts a first or¬ 
der transition from a dynamically stable symmetric so¬ 
lution to stable skewed solutions. In this calculation, T 
is taken to be 0.014—significantly higher values of this 
quantity lead to a washing out of evidence of a dynam¬ 
ical transition. Here the compliances are defined as in 
Sec. IV with x(t) replaced by (x(t)), the average being 
taken with respect to the distribution function P( x, t), 
i.e. (x(t)) = J^° oo x'P(x',t)dx'. The plot in Fig. 28 is to 


FIG. 27. One period (the twentieth) of the probability distri¬ 
bution, P(x, t ) as given by the solution of the Fokker-Planck 
equation ( |7.4[ )-( 7.61, with A = 2.7, uj = 0.5, k = 2, Fo = 0 
and F = 0.2. Also shown are the two stable skewed solu¬ 
tions (solid red curves) and the unstable symmetric solution 
(dashed red curve). Lighter colors correspond to higher values 
of P(x, t). 


I/A 



FIG. 28. The total compliance It/A (solid black curve), the 
dissipative compliance, Id/A (long-dashed blue curve) and the 
storage compliance, I s /A (dashed red curve). Here k = 2, 
to = 0.5, F 0 = 0 and T = 0.014. 


be compared with Fig. 16 the plot of compliances in the 
absence of noise. In the presence of noise the transition 
from an elastic to a viscous response regime no longer 
displays a first order discontinuity. 


VIII. CONCLUDING REMARKS 

The investigations of Zocchi et al. m of the re¬ 
sponse of a folded protein to a sinusoidal driving force 
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point to a description in terms of a single collective coor¬ 
dinate. Moreover, the response appears to be quantifiable 
in terms of a viscoelastic formalism, with the additional 
feature of a sharp transition from dominantly elastic re¬ 
sponse at low amplitude drive to mainly viscous response 
at high drive amplitude. We find that a simple model ex¬ 
hibits behavior consistent with those observations. From 
such analysis we conclude that the combination of experi¬ 
mental findings and characterization in terms of a simple 
model by Zocchi et al. forms a promising basis for a 
comprehensive description of the gross mechanical and 
dynamical properties of a class of folded proteins. Fur¬ 
thermore, we strongly believe that it ought to be useful 
in further investigations, both experimental and theoret¬ 
ical, of those properties. 
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Appendix A: The nonlinear restoring force 

The form for the restoring force derives from the func¬ 
tion 

A (x,a) = - o^o (Al) 

7r x z + or 


Extensions of the present study should aid in the de¬ 
velopment of insights into the structural properties of 
proteins, both in isolation and as components of larger 
systems. For example, it would be interesting to expand 
the characterization of a protein configuration beyond a 
single collective variable, to take into account internal 
adjustments to external influences. Given the size and 
complexity of those molecules the number of mechanical 
degrees of freedom will be considerable. However, ju¬ 
dicious analysis should allow the investigator to expand 
the set of dynamical variables to a manageable size. Ad¬ 
ditionally, the effects of inertia, while small, are almost 
certainly non-negligible, and the effects of noise deserve 
fuller attention. 


The existence of a dynamical phase transition in the 
simple model studied here deserves note in and of itself, 
and the genesis of this transition merits further study. 
It would be of great interest to identify the essential 
characteristics of a nonlinear system that lead to this 
phenomenon. As regards symmetry breaking, a crucial 
aspect of the model is the existence of a symmetry in 
the underlying equation of motion, in particular the fact 
that the restoring force, F(x ) satisfies F(—x) = —F(x); 
a similar symmetry is present in the Suzuki-Kubo model 
03132 ]. However, the protein conformation is almost cer¬ 
tainly not consistent with such a force, given that the 
collective co-ordinate x will have one sign—say positive— 
when the protein is stretched and the opposite sign when 
it is compressed. It is entirely reasonable to expect that 
the “cracking” process, in which a linear restoring force is 
replaced by another relationship between that force and 
the relative displacements in the protein, will differ in 
those two regimes, which means that there is no symme¬ 
try to break. Nevertheless, as shown in Appendix[F] even 
in the absence of such intrinsic asymmetry, the equation 
of motion can nevertheless lead to the kind of dynamical 
transition discussed in this paper. Thus, there are good 
reasons to expect that the response observed by Zocchi 
and co-workers does indeed arise from a true dynamical 
phase transition. 


In the limit a = 0, this is just the Dirac delta function. 
Two functions that result from integrating this function 
once and then twice are 


fi(x, a) 
h{x,a) 


2 tan 1 ( x/a ) 
7r 


(A2) 


2(xtan 1 {x/a) — \a In (a 2 + x 2 )) 


7T 

(A3) 


Then, with coefficients 

^ 2aFo — 7rfc(l + a 2 ) 

4 ((1 + a 2 ) tan -1 (1/a) — a) 
_Fo 
2 

2aFo — irk(l + a 2 ) 

4 ((1 + a 2 ) tan -1 (1/a) — a) 

The function 


(A4) 

(A5) 


A(/i(x + l,a) + fi(x - 1, a)) 

+B (f 2 (x - l,a) - f 2 (x+ l,a)) (A6) 

has the form shown in Fig. [2j with slope — k at the origin 
and asymptotes of ±Fq. The parameter a determines the 
sharpness of the transition from a linear restoring force to 
a restoring force that is constant. This cumbersome ex¬ 
pression for the restoring force enables independent vari¬ 
ation of the features controlled by a and Fq. Figure [29] 
displays three instances of the kind of restoring force that 
we explore. 

In the body of the paper we have taken the transition 
from Hookean behavior to occur around |rc| = 1. 


Appendix B: Pairs of skewed solutions 


The form of the equation of motion, (2.1) and the sym¬ 
metry of the restoring force, 


F(x) = —F(—x) 


(Bl) 
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F(x) 



FIG. 29. Three instances of restoring forces arising from (A 61 . 
In all cases, the absolute value of the slope of the curve at 
x = 0 is k = 2, and the absolute value of the asymptote is 
Fo = 0.3. The parameter a controlling the transition between 
the central and outer region is 0.01 for the solid curve, 0.1 
for the short-dashed curve and 0.5 for the long-dashed curve. 
The restoring force changes from Hookean to constant near 
x = ±1. 


or 


Sx(t + 2ir/(jj) 
5x{t) 


= exp 


nt-\-2ir/uj 


7 Jt 


F'(xo(t))dt 


(C3) 


Figure[30|shows a restoring force that changes abruptly 
from Hookean to vanishing, along with its derivative. In 


F(x), dF(x)/dx 



allow one to generate, from a solution x(t), a “mirror 
image” solution equa l to —x(t+ 7r/w). To see this, replace 
t by t + tt/u> in (2.1), which leads to 


dx(t + 7 t/lo) 
dt 


— [F(x(t + 7r/w)) + A sin(w(f + 7r/w))] 

7 

— [— F(—x(t + 7r/w)) — Asin(wf)] (B2) 
7 


If we now multiply both sides of the resulting equation by 
— 1, we end up with the result that —x{t + tt/uj) satisfies 
the same equation as x(t). This allows one to generate, 
from any skewed solution, a second one, skewed in the 
opposite direction. By contrast, a symmetric solution 
will reproduce itself under this transformation. 


Appendix C: Stability analysis for solutions to the 
equation of motion in Sec. m 


To assess the stability of a solution of the equation of 
motion when the restoring force exhibits a sharp break 
between Hooke’s Law and vanishing amplitude as in 
(3.1), we perturb that equation as follows. Assume a so¬ 


lution, Xq (t), to (2.1). Then, write x{t) = Xq( t) + Sx(t). 


At first order in the perturbation Sx(t), the equation be¬ 
comes 


^l = ±[6x(t)F'(x 0 m 


(Cl) 


Dividing both sides of (Cl) by Sx(t) and integrating over 
a period of the drive, we find 


In 


5x(t + 2 -k/oj) 
Sx(t) 


-f 

7 Jt 




F'(x 0 (t))dt (C2) 


FIG. 30. Restoring force with a rapid transition between 
Hooke’s law in a region about the origin (dashed curve) along 
with its derivative (solid curve). When the transition is per¬ 
fectly abrupt the peaks in the force derivative become Dirac 
delta functions. 


the limit of a discontinuous transition, the peaks in the 
derivative become Dirac delta functions. If we are inter¬ 
ested in the skewed solutions descri bed i n Sec. JHD then 
the integral on the right hand side of (C2) and (C3) yields 


three contributions. The first is from the interior region 
and will be equal to — k multiplied by the amount of time 
that x(t) < xq, which, according to (3.4), is 27 t/o 7 — 2t 0 - 


Then there is the result of integrating through the very 
narrow peaks in the derivative. Here, we write 


F'(xo{t))dt = 


F'(x 0 (t)) 


dxo(t)/dt 
F\x) 


dx 0 (t ) 


-> ± In 


F(x) ± Asin(o;to) 
Asin(o;to) 


dx 


Asin(wfo) T &%o 


(C4) 


The upper sign on the right hand side of (C4) applies 


in the case that the trajectory through the peak in the 
derivative is out of the inner region and the lower sign 
applies when the trajectory is into the inner region. The 
time, t has been set equal to to, as passage through the 
peaks occurs in effectively a single instant of time. As¬ 
sembling the two contributions we fnd for the right hand 
side of (C31 


exp 


k 

(2n \ 

1 , / 

— 

- 2t 0 

+ - In 

L 7 

K w / 

7 \ 


(.A sin(wfo)) 2 


(A sin(wfo)) 2 — (kx o) 2 


(C5) 
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If the expression in square brackets is real and negative 
then perturbations to the trajectory Xq (t) die off, and 
the periodic solution to the equation of motion is sta¬ 
ble. If it is positive, then those perturbations grow expo¬ 
nentially, and the solution is unstable. If the expression 
is imaginary, then the assumption of a solution is incor¬ 
rect, in that the velocity at one of the transitions between 
Hooke’s law force and no restoring force has the wrong 
sign. We find that in the case of graphical solutions to Eq. 
(3.8) that are not equal to 7 r/w, the one with a smaller 
value of to does not correspond to an actual trajectory, 
and the one with the larger value of to corresponds to a 
dynamically stable solution to the equation of motion. 


Appendix D: Review of the statistical mechanics of 
tricriticality 


minimum of the free energy is at ip = 0. On the other 
hand, when t < 0, there are two minima, spaced equidis¬ 
tant from ip = 0. The minimum at ip = 0 has become 
a free energy maximum. In mechanics, that maximum 
would correspond to a point of unstable equilibrium, in 
contrast to the minima, which represent points of stable 
equilibrium. However, a thermodynamic system will in¬ 
evitably fluctuate away from such a state. Because the 
free energy is even in ip, in that J-(—ip) = T{ip), each 
of the minima for t < 0, which lies to one side or the 
other of the origin, represents a violation of the symmetry 
of the physics underlying the system’s thermodynamics, 
and the transition that occurs as t passes through zero is 
called a symmetry breaking phase transition. The deter¬ 
mination of the minima follows from the solution of the 
equation 


According to the canonical model of tricritical points 
based on an order parameter, the thermodynamic behav¬ 
ior of a system near a symmetry breaking phase transi¬ 
tion is governed by a free energy of the form 

F{ip) =To + Ctip 2 + wip 4 + vip 6 (Dl) 


dT{ip) 

dip 


(D2) 


On the other hand, when u < 0, the behavior of T{ip) 
is shown in Fig. |32| In this case, the minima at non- 


Here, ip is the order parameter and t is a reduced temper¬ 
ature, proportional to the difference between the actual 
temperature of the system and a critical temperature. 
The coefficients C > 0, u and v are assumed to be insen¬ 
sitive to the temperature as long as t is sufficiently small. 
Thermodynamic stability requires v > 0; otherwise, the 
free energy takes on unbounded negative values as \ip\ 
increases. The system described by this free energy set¬ 
tles into the state of lowest free energy, determined by 
minimizing J-(ip). The results of this minimizing proce¬ 
dure depend qualitatively on the sign of the coefficient u. 
When u is positive, the ^-dependence of T is as shown in 
Fig. [3l] As shown in the figure, when t > 0, the unique 


Hi>) 



FIG. 31. Behavior of the free energy T(ip) in (Dll when the 
coefficients u and v are positive. The blue, dashed curve is 
the free energy when t > 0. The red, long-dashed, curve is 
the free energy when t < 0. The solid black curve is -F(ip) at 
the transition point, when t = 0. 





Ip 


FIG. 32. The graphical representation of cLF(ip)/dip. The 
curves follow the format of Fig. [31] 


zero values of ip develop while there is still a free energy 
minimum at ip = 0. As that occurs two free energy max¬ 
ima appear between the new minima and the free energy 
minimum at ip = 0. Eventually, at low enough temper¬ 
atures, the only minima are at finite ip, with one free 
energy maximum at ip = 0. Free energies and the related 
curves for the free energy derivative are displayed in Figs. 
[33] and [34] Focusing on the curves in Fig. [34] we see 
that at high enough temperatures, there is only one value 
of ip at which dF(ip)/dip passes through zero: ip = 0. As 
the temperature is lowered additional zeros appear, first 
four, two on each side of the origin, corresponding to a 
maximum and then a minimum in T(ip). At even lower 
temperatures, only three zeros remain, corresponding to 
minima at non-zero values of ip, while the zero at ip = 0 
refers to a free energy maximum at that point. 
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FIG. 33. Curves for J-(ip) when the fourth order coefficient u 
is negative. At high enough temperatures (blue long dashed 
curve) there is one free energy minimum, at ip = 0. As 
the temperature is lowered, two other minima emerge (blue 
dashed curve and solid black curve). In addition, two maxima 
that separate the minimum at ip = 0 and the flanking min¬ 
ima appear. At low enough temperatures the two minima at 
non-zero ip remain, the central extremum in the free energy 
having become a maximum (red long-dashed curve). 


dJ 7 (ip) / dip 



FIG. 34. Curves for dT(ip)/dip corresponding to the 
curves shown in Fig. [31] These curves pass through zero 
at extrema—maxima and minima—of J-(ip). They are for¬ 
matted in the same way as the free energy curves in Fig. |33| 


Appendix E: Alternate stability analysis of steady 
state solutions to <H3 


As an alternate approach to the stability analysis of 
steady state solutions to (6.11, obtai ned b y determining 
simultaneous solutions to (6.6) and (6.7), we return to 
the original equation (6.1) and expand about periodic 
solution, Xq( t). That is, write x(t) = Xo(t) +5x(t). The 
equation satisfied by Sx(t) is 


d 2 Sx(t) 

dt 2 


dSx(t) 


F' (x 0 (t))Sx(t) 


(El) 


Then, we write Sx(t) = ip(t)e 7t / 2 . The equation for 
ip(t) is 

- (l) 2 m = F'( Xo (t))m (E2) 

Because of the periodicity of Xq (t), this is just like the 
Schrodinger equation of a particle in a periodic potential. 
The periodicity is that of the driving force. We can now 
treat this as a problem in one dimensional band theory 
m We start with a solution of the above equation, 
ipi(t) that has the initial conditions 


V’i(O) = 1 (E3) 

V'i(O) = 0 (E4) 

and a function ip 2 (t) satisfying 

V> 2 (0) = 0 (E5) 

V4(0) = 1 (E6) 


A general solution will be a linear combination of ip\ (t) 
and ip 2 (t), where the coefficients are determined, for ex¬ 
ample, by the value of the solution and the value of its 
derivative at t = 0. If we write 


ip(t) = A'0i(t) + Bip 2 (t) 


(E7) 


Then, 

A = ip( 0) (E8) 

B = ip'( 0) (E9) 

After one full period, the magnitude and derivative of the 
solution is given by 

A' = Aipi(2ir / u>) + Bip 2 (2 t:/uj) (E10) 

B' = Aip[(2n/oj) + BiP[(2tt/uj) (Ell) 

This gives rise to the map 


A' 

B' 


/'0i(2tt/w) V , i(27t/w) 

V ^l( 27r / w ) V’2( 27r / W ) 


(E12) 


We know from the invariance of the Wronskian and from 
its value at t = 0 that 


i’idWid) - = 1 


(E13) 


for all t. This means that the determinant of the matrix 


in (E12) is equal to one. 


Given that the relationship (E12) holds for the values 


of A and B for each succeeding period, the ultimate be¬ 
havior of the solution is going to be determined by the 
eigenvalues of the matrix on the right hand side of that 
equation. In light of the relationship between 5x(t) and 
ip(t ), we see that the stability of the solution is going to 
be determined by the exponential e _77r// “ and the eigen¬ 
values. Specifically, by the product 


0 —77r/u 


iPi(2tt/uj) + ip' 2 (2'K/ui) 


± 


\j W>i(2t r/w) + Ip’ 2 (2-K/u)) 2 - 4 


(E14) 


dt 
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If the absolute value of (E14) is less than one, for both 
signs in the expression, then the solution is stable, and if 
that is not the case the solution is unstable. The result in 
(E14) follows in part from the fact that the determinant 
of the matrix in (E12) is equal to one. 

Making use of the above method, we arrive at the same 
conclusions regarding the steady state solutions to the 


equation of motion (6.1) shown in Sec. VI 


Appendix F: Dynamical transition when the force is 
inherently asymmetric 


We have been considering a system in which the force, 
F(x) that restores it to stable equilibrium in the absence 
of an external drive is symmetric in the dynamical vari¬ 
able that characterizes its distortion from that state in 
that F(—x) = —F(x). It is natural to ask whether this 
symmetry in the equation of motion is essential to the 
existence of the dynamical transitions discussed here. As 
it turns out those transitions can also arise when the un¬ 
derlying symmetry is absent. Consider the restoring force 
shown in Fig. [35] As shown in that figure, the force is 


F(x) 



FIG. 35. A restoring force, F(x) that is Hookean for all posi¬ 
tive values of x but that reverts to zero at x = — 1 . The spring 
constant in the Hookean regime is k = 2. 


x(2jr/co) 
x(0) 



FIG. 36. The map from x(0) to x(2n/u ;) generated by the 
equation of motion (2.1l with the restoring force shown in 
Fig. [35] and u> = 2, 7 = 1 and A = 2.762. The points of 
intersection between the map and the 45° line, also shown, 
are indicated by small circles. The two outer intersections 
correspond to dynamically stable steady state solutions to the 
equation of motion, while the inner intersection corresponds 
to a dynamically unstable steady state solution. 


x(t) 



of the form F = —kx, with spring constant k = 2 for 
x < 1. Outside of that region the restoring force is equal 
to zero. 

In Fig. 36 we see the map from x(0) to x(2tt/uj), along 
with the 45° line for drive amplitude A = 2.762 , fre - 
quency ui = 2 and 7 = 1 in the equation of motion ( 2 . 11 . 
This kind of map is consistent with the existence of two 
stable fixed point solutions and one unstable solution; 
Fig, [37] shows those three solutions. This coexistence 
points to a dynamical transition identical to sharp tran¬ 
sitions described in the body of this paper, the only dif¬ 
ference being the absence of symmetry breaking; there is 
no underlying symmetry to break. 

As further evidence for the dynamical transition, Fig. 
[38] contains the maps for four different drive amplitudes 
in the same region plotted in Fig. [37| This plot is to 


FIG. 37. The three steady state solutions to the equation of 
motion corresponding to the points of intersection in Fig. |36| 
plotted over two periods of the drive, Asin(u;t). The thick 
red curve corresponds to a trajectory that lies entirely in the 
region in which the restoring force is linear. The thin blue 
curve is the stable solution that extends outside that region 
(,x < —1), and the dashed curve is the dynamically unstable 
solution that separates the two stable solutions for x(t). 


be compared with the maps shown in Fig. |17| The pa¬ 
rameters are as described in Fig. [36j except for the drive 
amplitudes, which range from A = 2.755 to A = 2.78. 
The discussion of the scenario corresponding to this fig¬ 
ure parallels the corresponding commentary just above 

Fig. [17] 
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x(2nla>) 



The noteworthy results here are, first, that at the high¬ 
est drive amplitude there is only one skewed steady state 
solution to the equation of motion instead of two and, 
second, that when there is more than one solution, the 
two dynamically stable ones correspond to a symmetric 
solution and a single skewed solution, separated by an 
unstable steady state solution, as plotted in Fig. [37} 


FIG. 38. Illustrating a first order transition in the case in 
which the restoring force is as shown in Fig. |35[ T he pa¬ 
rameters are as described in the caption to Fig. |36| except 
for the values of the drive amplitudes. Those amplitudes are 
2.755 (intermediate dashed blue curve), 2.762 (short dashed 
red curve), 2.775 (solid black curve) and 2.78 (long dashed red 
curve). Note that there are no more than three steady state 
solutions to the equation of motion and that at the extremes 
there is only one steady state solution to the equation of mo¬ 
tion o in contrast to the case of the symmetry breaking 
dynamical phase transition. 


[1] Y. Wang and G. Zocchi, Plos One 6 (2011). 

[2] Y. Wang and G. Zocchi, EPL 96 (2011). 

[3] W. Yong and G. Zocchi, Europhysics Letters 96, 18003 

(6 PP .) ( 2011 ). 

[4] Q. Hao, J. Landy, and G. Zocchi, Physical Review 
E (Statistical, Nonlinear, and Soft Matter Physics) 86, 
041915 (2012). 

[5] D. Whitford, Proteins : structure and function (J. Wiley 
& Sons, Hoboken, NJ, 2005) david Whitford. 

[6] R. Phillips, J. Kondev, J. Theriot, and H. G. Garcia., 
Physical biology of the cell, 2nd ed. (Garland Science, 
London New York, NY, 2013). 

[7] R. S. Lakes, Viscoelastic materials (Cambridge Univer¬ 
sity Press, Cambridge; New York, 2009) pp. xviii, 461 
P- 

[8] M. Suzuki and R. Kubo, Journal of the Physical Society 
of Japan 24 , 51 (1968) 

[9] T. Tome and M. J. de Oliveira, Physical Review A 41 , 
4251 (1990) 

[10] W. N. Findley, J. S. Lai, and K. Onaran, Creep and 
relaxation of nonlinear viscoelastic materials : with an 
introduction to linear viscoelasticity, Dover books on en¬ 
gineering (Dover, New York, 1989) pp. xii, 371 p. 


[11] M. M. Tirion, Physical Review Letters 77 , 1905 (1996). 

[12] A. Ansari, J. Berendzen, S. F. Bowne, H. Frauenfclder, 
I. E. T. Iben, T. B. Sauke, E. Shyamsunder, and R. D. 
Young, Proceedings of the National Academy of Sciences 
of the United States of America 82, 5000 (1985) 

[13] O. Miyashita, J. N. Onuchic, and P. G. Wolynes, Pro¬ 
ceedings of the National Academy of Sciences of the 
United States of America 100, 12570 (2003). 

[14] See also Anderson et. al: Chapter 8 of A. Noy, Handbook 
of molecular force spectroscopy (Springer, New York, NY, 
2008). 

[15] J. Guckenheimer and P. Holmes, Nonlinear oscillations, 
dynamical systems, and bifurcations of vector fields, corr. 
5th print, ed., Applied mathematical sciences (Springer, 
New York, 1997). 

[16] M. J. Feigenbaum, Journal of Statistical Physics 19 , 25 
(1978) 

[17] H. Risken, The Fokker-Planck equation : methods of so¬ 
lution and applications, 2nd ed., Springer series in syner¬ 
getics, (Springer-Verlag, Berlin ; New York, 1996). 

[18] N. W. Ashcroft and N. D. Mermin, Solid state physics 
(Holt Rinehart and Winston, New York,, 1976) pp. xxi, 
826. 











